library(lme4)
## Loading required package: Matrix
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(stringr)
library("wesanderson")
library(ggnewscale)
library('corrr')
library('ggcorrplot')
library("FactoMineR")
library("factoextra")
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
setwd("~/Downloads/QueenModel_052624")
prefixes = c("RooibosTea_QR_1216_1646", "RooibosTea_QL_1216_1646", "MexHotChoc_QR_1216_1646", "MexHotChoc_QL_1216_1646", "20230213_1745_AlmdudlerGspritzt_C1", "20230213_1745_AlmdudlerGspritzt_C0", "20221209_1613_QR", "20221209_1613_QL", "20221123_1543_AmericanoLatte_QR", "20221123_1543_AmericanoLatte_QL")
Day = 1
Days <- list()
Start = 0
#Total Degree
for(i in 1:10){
for(j in 0:97){
if(file.exists(paste(prefixes[i],sprintf("%03d", j),"Head_to_Head_False_Centralities.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],sprintf("%03d", j),"Head_to_Head_False_Centralities.csv",sep = "_"))
Deg$Node = sub('\\.0','',Deg$X)
Deg$Day = floor((j/24)+1)
Deg$Hour = sprintf("%03d", j)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$ID = paste(Deg$Col, Deg$Node, Deg$Hour, sep='_')
Deg$Bee = paste(Deg$Col, Deg$Node, sep='_')
Deg$Mod = "Head_Head"
Deg$Trial = str_extract(Deg$Col, ".+?(?=_)")
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalDeg = rbind(TotalDeg, Deg)
}
if(Start == 0){
TotalDeg = Deg
Start = 1
}
}
}
}
names(TotalDeg)[names(TotalDeg) == "degree"] <- "Degree"
names(TotalDeg)[names(TotalDeg) == "between"] <- "Between"
names(TotalDeg)[names(TotalDeg) == "close"] <- "Close"
names(TotalDeg)[names(TotalDeg) == "eigen"] <- "Eigen"
#Bout Degree
Start = 0
for(i in 1:10){
for(j in 0:97){
if(file.exists(paste(prefixes[i],sprintf("%03d", j),"Head_to_Head_True_Centralities.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],sprintf("%03d", j),"Head_to_Head_True_Centralities.csv",sep = "_"))
Deg$Node = sub('\\.0','',Deg$X)
Deg$Day = floor((j/24)+1)
Deg$Hour = sprintf("%03d", j)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$ID = paste(Deg$Col, Deg$Node, Deg$Hour, sep='_')
Deg$Bee = paste(Deg$Col, Deg$Node, sep='_')
Deg$Mod = "Head_Head"
Deg$Trial = str_extract(Deg$Col, ".+?(?=_)")
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalBout = rbind(TotalBout, Deg)
}
if(Start == 0){
TotalBout = Deg
Start = 1
}
}
}
}
TotalBoutSubset <- TotalBout[, c("ID", "degree", "between", "close", "eigen")]
names(TotalBoutSubset)[names(TotalBoutSubset) == "degree"] <- "boutDegree"
names(TotalBoutSubset)[names(TotalBoutSubset) == "between"] <- "boutBetween"
names(TotalBoutSubset)[names(TotalBoutSubset) == "close"] <- "boutClose"
names(TotalBoutSubset)[names(TotalBoutSubset) == "eigen"] <- "boutEigen"
BigDataSheet <- merge(TotalDeg, TotalBoutSubset, by = "ID")
#Body Degree
Start = 0
for(i in 1:10){
for(j in 0:97){
if(file.exists(paste(prefixes[i],sprintf("%03d", j),"Head_to_Body_False_Centralities.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],sprintf("%03d", j),"Head_to_Body_False_Centralities.csv",sep = "_"))
Deg$Node = sub('\\.0','',Deg$X)
Deg$Day = floor((j/24)+1)
Deg$Hour = sprintf("%03d", j)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$ID = paste(Deg$Col, Deg$Node, Deg$Hour, sep='_')
Deg$Bee = paste(Deg$Col, Deg$Node, sep='_')
Deg$Mod = "Head_Head"
Deg$Trial = str_extract(Deg$Col, ".+?(?=_)")
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalBody = rbind(TotalBody, Deg)
}
if(Start == 0){
TotalBody = Deg
Start = 1
}
}
}
}
TotalBodySubset <- TotalBody[, c("ID", "degree", "between", "close", "eigen")]
names(TotalBodySubset)[names(TotalBodySubset) == "degree"] <- "bodyDegree"
names(TotalBodySubset)[names(TotalBodySubset) == "between"] <- "bodyBetween"
names(TotalBodySubset)[names(TotalBodySubset) == "close"] <- "bodyClose"
names(TotalBodySubset)[names(TotalBodySubset) == "eigen"] <- "bodyEigen"
BigDataSheet <- merge(BigDataSheet, TotalBodySubset, by = "ID")
#Bout Length
BigDataSheet$AverageBoutLength = BigDataSheet$Degree / BigDataSheet$bodyDegree
#Ovary Measurements
Ovaries = read.csv('OvaryMeasurements.csv')
Ovaries$AverageOvaryLength = (Ovaries$LongestOocyteLength1..mm. + Ovaries$LongestOocyteLength2..mm.) / 2
Ovaries$AverageOvaryWidth = (Ovaries$LongestOocyteWidth1..mm. + Ovaries$LongestOocyteWidth2..mm.) / 2
Ovaries$Bee = paste(Ovaries$ColonyID, '_ArUcoTag#', Ovaries$Tag, sep='')
BigDataSheet <- merge(BigDataSheet, Ovaries[, c("Bee", "AverageOvaryLength", "AverageOvaryWidth")], by = "Bee", all.x = TRUE)
#Wing Measurements
Wings = read.csv('WingMeasurements.csv')
Wings$AverageWingLength = (Wings$wing_1_mm + Wings$wing_2_mm) / 2
Wings$TagNumber = str_extract(Wings$ID, "(?<=Tag)\\d+")
Wings$Bee = paste(Wings$ColonyID, '_ArUcoTag#', Wings$TagNumber, sep='')
BigDataSheet <- merge(BigDataSheet, Wings[, c("Bee", "AverageWingLength")], by = "Bee", all.x = TRUE)
#Track Presence
Start=0
for(i in 1:10){
for(j in 0:97){
if(file.exists(paste(prefixes[i],sprintf("%03d", j),"8_8_3_1_10_interp_filtered.h5_TrackPresence.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],sprintf("%03d", j),"8_8_3_1_10_interp_filtered.h5_TrackPresence.csv",sep = "_"))
Deg$track = sub('\\.0','',Deg$track)
Deg$Day = floor((j/24)+1)
Deg$Hour = sprintf("%03d", j)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$ID = paste(Deg$Col, Deg$track, Deg$Hour, sep='_')
Deg$Bee = paste(Deg$Col, Deg$track, sep='_')
Deg$Mod = "Head_Head"
Deg$Trial = str_extract(Deg$Col, ".+?(?=_)")
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalPresence = rbind(TotalPresence, Deg)
}
if(Start == 0){
TotalPresence = Deg
Start = 1
}
}
}
}
TotalPresenceSubset <- TotalPresence[, c("ID", "Presence")]
BigDataSheet <- merge(BigDataSheet, TotalPresenceSubset, by = "ID")
#Antennna Presence
Start=0
for(i in 1:10){
for(j in 0:97){
if(file.exists(paste(prefixes[i],sprintf("%03d", j),"8_8_3_1_10_interp_filtered.h5_TrackPresenceAnt.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],sprintf("%03d", j),"8_8_3_1_10_interp_filtered.h5_TrackPresenceAnt.csv",sep = "_"))
Deg$track = sub('\\.0','',Deg$track)
Deg$Day = floor((j/24)+1)
Deg$Hour = sprintf("%03d", j)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$ID = paste(Deg$Col, Deg$track, Deg$Hour, sep='_')
Deg$Bee = paste(Deg$Col, Deg$track, sep='_')
Deg$Mod = "Head_Head"
Deg$Trial = str_extract(Deg$Col, ".+?(?=_)")
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalPresence = rbind(TotalPresence, Deg)
}
if(Start == 0){
TotalPresence = Deg
Start = 1
}
}
}
}
TotalPresenceSubset <- TotalPresence[, c("ID", "Presence")]
names(TotalPresenceSubset)[names(TotalPresenceSubset) == "Presence"] <- "AntPresence"
BigDataSheet <- merge(BigDataSheet, TotalPresenceSubset, by = "ID")
#Vels Data
Start=0
for(i in 1:10){
for(j in 0:97){
if(file.exists(paste(prefixes[i],sprintf("%03d", j),"velstrackqueen.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],sprintf("%03d", j),"velstrackqueen.csv",sep = "_"))
Deg$X = sub('\\.0','',Deg$X)
Deg$Day = floor((j/24)+1)
Deg$Hour = sprintf("%03d", j)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$ID = paste(Deg$Col, Deg$X, Deg$Hour, sep='_')
Deg$Bee = paste(Deg$Col, Deg$X, sep='_')
Deg$Mod = "Head_Head"
Deg$Trial = str_extract(Deg$Col, ".+?(?=_)")
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalVels = rbind(TotalVels, Deg)
}
if(Start == 0){
TotalVels = Deg
Start = 1
}
}
}
}
TotalVelsSubset <- TotalVels[, c("ID", "mean_vel", "move_perc")]
BigDataSheet <- merge(BigDataSheet, TotalVelsSubset, by = "ID")
Start=0
for(i in 1:10){
if(file.exists(paste(prefixes[i],"096_DispersionMetrics.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],"096_DispersionMetrics.csv",sep = "_"))
Deg$Tag = sub('\\.0','',Deg$Tag)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$Bee = paste(Deg$Col, Deg$Tag, sep='_')
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalDisp = rbind(TotalDisp, Deg)
}
if(Start == 0){
TotalDisp = Deg
Start = 1
}
}
}
DispNew <- reshape(TotalDisp, idvar=c("Tag", "Col", "QR", "Bee", "Queen"), timevar = "Day", direction="wide")
DispSubset <- DispNew[, c("Bee", "N90.Day1", "N90.Day2", "N90.Day3", "N90.Day4", "MRSD.Day1", "MRSD.Day2", "MRSD.Day3", "MRSD.Day4")]
BigDataSheet <- merge(BigDataSheet, DispSubset, by = "Bee", all.x = TRUE)
#Diectional Degree
Start = 0
for(i in 1:10){
for(j in 0:97){
if(file.exists(paste(prefixes[i],sprintf("%03d", j),"DegDir.csv",sep = "_"))) {
Deg = read.csv(paste(prefixes[i],sprintf("%03d", j),"DegDir.csv",sep = "_"))
Deg$Node = sub('\\.0','',Deg$Origin.interactor)
Deg$Day = floor((j/24)+1)
Deg$Hour = sprintf("%03d", j)
Deg$Col = prefixes[i]
Deg$QR = i %% 2 == 1
Deg$ID = paste(Deg$Col, Deg$Node, Deg$Hour, sep='_')
Deg$Bee = paste(Deg$Col, Deg$Node, sep='_')
Deg$Mod = "Head_Head"
Deg$Trial = str_extract(Deg$Col, ".+?(?=_)")
Deg$Queen = Deg$Bee %in% c("RooibosTea_QR_1216_1646_ArUcoTag#52", "MexHotChoc_QR_1216_1646_ArUcoTag#13", "20230213_1745_AlmdudlerGspritzt_C1_ArUcoTag#24", "20221209_1613_QR_ArUcoTag#47", "20221123_1543_AmericanoLatte_QR_ArUcoTag#43")
if(Start == 1){
TotalDir = rbind(TotalDir, Deg)
}
if(Start == 0){
TotalDir = Deg
Start = 1
}
}
}
}
DirSubset <- TotalDir[, c("ID", "count_ori", "count_des", "Initiation.Freq")]
BigDataSheet <- merge(BigDataSheet, DirSubset, by = "ID")
BDSMeans <- aggregate(cbind(Degree,Close,Eigen,Between,QR,Queen,boutDegree,boutBetween,boutClose,boutEigen,bodyDegree,bodyBetween,bodyClose,bodyEigen,AverageBoutLength,AverageOvaryWidth,AverageWingLength,Presence,AntPresence,mean_vel,move_perc,N90.Day4,MRSD.Day4, Initiation.Freq) ~ Bee, BigDataSheet, mean)
BDSMeans$Trial = str_extract(BDSMeans$Bee, ".+?(?=_)")
BDSMeans$Col = str_extract(BDSMeans$Bee, ".+?(?=#)")
BDSMeansNoOv <- aggregate(cbind(Degree,Close,Eigen,Between,QR,Queen,boutDegree,boutBetween,boutClose,boutEigen,bodyDegree,bodyBetween,bodyClose,bodyEigen,AverageBoutLength,Presence,AntPresence,mean_vel,move_perc,N90.Day4,MRSD.Day4, Initiation.Freq) ~ Bee, BigDataSheet, mean)
BDSMeansNoOv$Trial = str_extract(BDSMeansNoOv$Bee, ".+?(?=_)")
BDSMeansNoOv <- BDSMeansNoOv %>%
mutate(QR_Queen_Condition = case_when(
QR==0 & Queen==0 ~ "Queenless",
QR==1 & Queen==0 ~ "Queenright",
Queen==1 ~ "Queen",
TRUE ~ NA_character_ # This handles any other case, which shouldn't exist in your scenario
)) %>%
mutate(QR_Queen_Condition = factor(QR_Queen_Condition, levels = c( "Queenright","Queenless", "Queen")))
BDSMeansNoOv$ID = paste(BDSMeansNoOv$Trial,BDSMeansNoOv$QR_Queen_Condition)
BDSMeanOfMean <- aggregate(cbind(Degree,Close,Eigen,Between,QR,Queen,boutDegree,boutBetween,boutClose,boutEigen,bodyDegree,bodyBetween,bodyClose,bodyEigen,AverageBoutLength,Presence,AntPresence,mean_vel,move_perc,N90.Day4,MRSD.Day4, Initiation.Freq) ~ ID, BDSMeansNoOv, mean)
BDSMeanOfMean$Trial = str_extract(BDSMeanOfMean$ID, ".+?(?= )")
BDSMeanOfMean <- BDSMeanOfMean %>%
mutate(QR_Queen_Condition = case_when(
QR==0 & Queen==0 ~ "Queenless",
QR==1 & Queen==0 ~ "Queenright",
Queen==1 ~ "Queen",
TRUE ~ NA_character_ # This handles any other case, which shouldn't exist in your scenario
)) %>%
mutate(QR_Queen_Condition = factor(QR_Queen_Condition, levels = c( "Queenright","Queenless", "Queen")))
ggplot(BDSMeanOfMean, aes(x = fct_rev(as.factor(QR_Queen_Condition)), y = Degree)) +
geom_line(aes(group = Trial), color = "darkgray") +
geom_point(aes(color = Trial), size = 5) +
scale_color_manual(values = wes_palette("Cavalcanti1")) +
xlab("") +
labs(color = "Source Colony") +
ylab("Mean Number of Interactions per Hour") + # Adjust axis labels
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "none",
text = element_text(size = 16),
panel.grid.major.x = element_line(color = "grey", linetype = "dashed"), # Keep vertical grid lines
panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"), # Keep vertical grid lines
panel.grid.major.y = element_blank(), # Remove horizontal grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal grid lines
axis.line.x = element_line(color = "black", size = 0.5), # Add x-axis line
axis.line.y = element_line(color = "black", size = 0.5), # Add y-axis line
strip.text = element_text(size = 14, face = "bold"), # Make facet plot titles larger and bold
axis.text.y.right = element_blank(), # Remove right y-axis text
axis.ticks.y.right = element_blank(), # Remove right y-axis ticks
aspect.ratio = 1
) +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5))
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
###MS:Figure 2b-----
BigDataSheetF2 <- BigDataSheet %>%
mutate(QR_Queen_Condition = case_when(
!QR & !Queen ~ "Queenless",
QR & !Queen ~ "Queenright",
Queen ~ "Queen",
TRUE ~ NA_character_ # This handles any other case, which shouldn't exist in your scenario
)) %>%
mutate(QR_Queen_Condition = factor(QR_Queen_Condition, levels = c( "Queenright","Queenless", "Queen")))
BigDataSheetF2$Alpha <- ifelse(BigDataSheetF2$Queen, .5, 0.005)
BigDataSheetF2 <- BigDataSheetF2[sample(nrow(BigDataSheetF2)), ]
BigDataSheetF2$PointSize <- ifelse(BigDataSheetF2$Queen, .003, .001)
ggplot(BigDataSheetF2, aes(x = as.integer(Hour), y = Degree/20, group = ID)) +
geom_jitter(aes(color = QR_Queen_Condition, alpha = Alpha, size = PointSize)) +
geom_smooth(aes(group = QR_Queen_Condition, color = QR_Queen_Condition)) +
scale_size(range = c(.001, .2)) +
scale_color_manual(
labels = c( "Queenright Worker", "Queenless Worker", "Queen"),
values = c( "#429CF0", "#161414","#7851A9"),
guide = guide_legend(direction = "horizontal")
) +
scale_x_continuous(breaks = c(0, seq(24, 96, by = 24)), limits = c(0, NA), expand = c(0, 0)) + # Expand limits to include 0
labs(color = "") +
xlab("Hour") +
ylab("Standardized Time Spent Interacting per Hour (s)") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = c(1, 1.02), # Move legend to top right
legend.justification = c(1, 1), # Align legend to top right
legend.background = element_rect(fill = alpha("white", 1)),
legend.title = element_text(size = 8, face = "bold"),
legend.text = element_text(size = 7),
panel.grid.major.x = element_line(color = "grey", linetype = "dashed"), # Keep vertical grid lines
panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"), # Keep vertical grid lines
panel.grid.major.y = element_blank(), # Remove horizontal grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal grid lines
axis.line.x = element_line(color = "black", size = 0.5), # Add x-axis line
axis.line.y = element_line(color = "black", size = 0.5), # Add y-axis line
strip.text = element_text(size = 10, face = "bold"), # Make facet plot titles larger and bold
panel.grid = element_blank(),
panel.border = element_blank()
) +
# scale_y_log10() +
guides(alpha = "none", size = "none")
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## Warning: Removed 261 rows containing missing values or values outside the scale range
## (`geom_point()`).
###MS:Figure 3------------
numerical_data <- BigDataSheet[,c(4:7,16:42)]
numerical_data = numerical_data[complete.cases(numerical_data), ]
data_normalized <- scale(numerical_data)
data.pca <- princomp(data_normalized)
fviz_eig(data.pca, addlabels = TRUE)
fviz_pca_var(data.pca, col.var = "black")
numerical_data$QR = BigDataSheet[complete.cases(BigDataSheet), ]$QR
numerical_data$Trial = BigDataSheet[complete.cases(BigDataSheet), ]$Trial
fviz_pca_ind(data.pca, label="none", habillage=numerical_data$QR,
addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2")
fviz_pca_ind(data.pca, label="none", habillage=numerical_data$Trial,
addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2")
numerical_data <- BDSMeans[,c(2:5,8:25)]
numerical_data = numerical_data[complete.cases(numerical_data), ]
data_normalized <- scale(numerical_data)
data.pca <- princomp(data_normalized)
fviz_eig(data.pca, addlabels = TRUE)
fviz_pca_var(data.pca, col.var = "black")
numerical_data$QR = BDSMeans[complete.cases(BDSMeans), ]$QR
numerical_data$Trial = BDSMeans[complete.cases(BDSMeans), ]$Trial
numerical_data$Queen = BDSMeans[complete.cases(BDSMeans), ]$Queen
numerical_data$OvaryIndex = numerical_data$AverageOvaryWidth/numerical_data$AverageWingLength
fviz_pca_ind(data.pca, label="none", habillage=numerical_data$QR,
addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2") + theme_minimal()
fviz_pca_ind(data.pca, label="none", habillage=numerical_data$Trial,
addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2")
fviz_pca_ind(data.pca, label="none", habillage=numerical_data$Queen,
addEllipses=TRUE, ellipse.level=0.95, palette = "Dark2")
fviz_pca_ind(data.pca, col.ind = numerical_data$OvaryIndex,geom="point",
legend.title = "OvaryIndex",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
BDSMeansQR = BDSMeans[BDSMeans$QR == TRUE,]
BDSMeansQL = BDSMeans[BDSMeans$QR == FALSE,]
cor.test(BDSMeansQR$Between, BDSMeansQR$boutBetween)
##
## Pearson's product-moment correlation
##
## data: BDSMeansQR$Between and BDSMeansQR$boutBetween
## t = 12.511, df = 126, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6554497 0.8128603
## sample estimates:
## cor
## 0.7443208
cor.test(BDSMeansQL$Between, BDSMeansQL$boutBetween)
##
## Pearson's product-moment correlation
##
## data: BDSMeansQL$Between and BDSMeansQL$boutBetween
## t = 28.646, df = 127, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9029735 0.9505362
## sample estimates:
## cor
## 0.9305802
###MS:Figure 4-------
BDSMeans$Alpha <- ifelse(BDSMeans$Queen, 0.5, 0.75)
BDSMeans$PointSize <- ifelse(BDSMeans$Queen, 2, 1)
ggplot(BDSMeans, aes(x = move_perc, y = Degree, group = interaction(QR, Queen))) +
geom_smooth(
data = subset(BDSMeans, Queen == FALSE), # Subset the data to exclude queens
method = lm, # Change method to lm for linear model
se = T, # Standard error
size = .75, # Adjust size to match first plot
linetype = 1, # Line type
aes(color = interaction(QR, Queen))
) + # Color by treatment
new_scale_color() +
geom_point(aes(color = Trial,shape=as.factor(QR)), alpha = 1,size = BDSMeans$PointSize) +
scale_color_manual(values = wes_palette("Cavalcanti1")) +
scale_x_continuous() + # Use a continuous scale for x
scale_y_continuous() + # Use a continuous scale for y
# labs(title="Total Interaction Count vs Movement Percentage", color = "Treatment") + # Adjust title and legend title
xlab("Fraction of Time Spent Moving") +
ylab("Mean Number of Interactions per Hour") + # Adjust axis labels
theme_minimal() +
theme(
text = element_text(size = 16),
plot.title = element_text(hjust = 0.5),
legend.position = "none",
panel.grid.major.x = element_line(color = "grey", linetype = "dashed"), # Keep vertical grid lines
panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"), # Keep vertical grid lines
panel.grid.major.y = element_blank(), # Remove horizontal grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal grid lines
axis.line.x = element_line(color = "black", size = 0.5), # Add x-axis line
axis.line.y = element_line(color = "black", size = 0.5), # Add y-axis line
strip.text = element_text(size = 14, face = "bold"), # Make facet plot titles larger and bold
axis.text.y.right = element_blank(), # Remove right y-axis text
axis.ticks.y.right = element_blank(), # Remove right y-axis ticks
aspect.ratio = 1
) +
xlim(0,1)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using formula = 'y ~ x'
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## `geom_smooth()` using formula = 'y ~ x'
BDSMeans$Alpha <- ifelse(BDSMeans$Queen, 0.5, 0.75)
BDSMeans$PointSize <- ifelse(BDSMeans$Queen, 2, 1)
BDSMeans$OvaryIndex = BDSMeans$AverageOvaryWidth/BDSMeans$AverageWingLength
ggplot(BDSMeans, aes(x = OvaryIndex, y = Degree, group = interaction(QR, Queen))) +
geom_smooth(
data = subset(BDSMeans, Queen == FALSE), # Subset the data to exclude queens
method = lm, # Change method to lm for linear model
se = T, # Standard error
size = .75, # Adjust size to match first plot
linetype = 1, # Line type
aes(color = interaction(QR, Queen))
) + # Color by treatment
new_scale_color() +
geom_point(aes(color = Trial,shape=as.factor(QR)), alpha = 1,size = BDSMeans$PointSize) +
scale_color_manual(values = wes_palette("Cavalcanti1")) +
scale_x_continuous() + # Use a continuous scale for x
scale_y_continuous() + # Use a continuous scale for y
# labs(title="Total Interaction Count vs Movement Percentage", color = "Treatment") + # Adjust title and legend title
xlab("Ovary Index") +
ylab("Mean Number of Interactions per Hour") + # Adjust axis labels
theme_minimal() +
theme(
text = element_text(size = 16),
plot.title = element_text(hjust = 0.5),
legend.position = "none",
panel.grid.major.x = element_line(color = "grey", linetype = "dashed"), # Keep vertical grid lines
panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"), # Keep vertical grid lines
panel.grid.major.y = element_blank(), # Remove horizontal grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal grid lines
axis.line.x = element_line(color = "black", size = 0.5), # Add x-axis line
axis.line.y = element_line(color = "black", size = 0.5), # Add y-axis line
strip.text = element_text(size = 14, face = "bold"), # Make facet plot titles larger and bold
axis.text.y.right = element_blank(), # Remove right y-axis text
axis.ticks.y.right = element_blank(), # Remove right y-axis ticks
aspect.ratio = 1
) +
xlim(0,1)
## `geom_smooth()` using formula = 'y ~ x'
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## `geom_smooth()` using formula = 'y ~ x'
BigDataSheetF2$StandardizedDegree = BigDataSheetF2$Degree * BigDataSheetF2$AntPresence
ggplot(BigDataSheetF2, aes(x = as.integer(Hour), y = StandardizedDegree/20, group = ID)) +
geom_jitter(aes(color = QR_Queen_Condition, alpha = Alpha, size = PointSize)) +
geom_smooth(aes(group = QR_Queen_Condition, color = QR_Queen_Condition)) +
scale_size(range = c(.001, .2)) +
scale_color_manual(
labels = c( "Queenright Worker", "Queenless Worker", "Queen"),
values = c( "#429CF0", "#161414","#7851A9"),
guide = guide_legend(direction = "horizontal")
) +
scale_x_continuous(breaks = c(0, seq(24, 96, by = 24)), limits = c(0, NA), expand = c(0, 0)) + # Expand limits to include 0
labs(color = "") +
xlab("Hour") +
ylab("Standardized Time Spent Interacting per Hour (s)") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = c(1, 1.02), # Move legend to top right
legend.justification = c(1, 1), # Align legend to top right
legend.background = element_rect(fill = alpha("white", 1)),
legend.title = element_text(size = 8, face = "bold"),
legend.text = element_text(size = 7),
panel.grid.major.x = element_line(color = "grey", linetype = "dashed"), # Keep vertical grid lines
panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"), # Keep vertical grid lines
panel.grid.major.y = element_blank(), # Remove horizontal grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal grid lines
axis.line.x = element_line(color = "black", size = 0.5), # Add x-axis line
axis.line.y = element_line(color = "black", size = 0.5), # Add y-axis line
strip.text = element_text(size = 10, face = "bold"), # Make facet plot titles larger and bold
panel.grid = element_blank(),
panel.border = element_blank()
) +
# scale_y_log10() +
guides(alpha = "none", size = "none")
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
## Warning: Removed 242 rows containing missing values or values outside the scale range
## (`geom_point()`).
###MS:Supp Figure 5--------
ggplot(BDSMeanOfMean, aes(x = fct_rev(as.factor(QR_Queen_Condition)), y = N90.Day4)) +
geom_line(aes(group = Trial), color = "darkgray") +
geom_point(aes(color = Trial), size = 5) +
scale_color_manual(values = wes_palette("Cavalcanti1")) +
xlab("") +
labs(color = "Source Colony") +
ylab("N90 (Dispersion Measure)") + # Adjust axis labels
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "none",
text = element_text(size = 16),
panel.grid.major.x = element_line(color = "grey", linetype = "dashed"), # Keep vertical grid lines
panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"), # Keep vertical grid lines
panel.grid.major.y = element_blank(), # Remove horizontal grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal grid lines
axis.line.x = element_line(color = "black", size = 0.5), # Add x-axis line
axis.line.y = element_line(color = "black", size = 0.5), # Add y-axis line
strip.text = element_text(size = 14, face = "bold"), # Make facet plot titles larger and bold
axis.text.y.right = element_blank(), # Remove right y-axis text
axis.ticks.y.right = element_blank(), # Remove right y-axis ticks
aspect.ratio = 1
) +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5))
###Presence Distribution-----
BDSRanked <- aggregate(cbind(Degree,Close,Eigen,Between,QR,Queen,boutDegree,boutBetween,boutClose,boutEigen,bodyDegree,bodyBetween,bodyClose,bodyEigen,AverageBoutLength,Presence,AntPresence,mean_vel,move_perc,N90.Day4,MRSD.Day4,Initiation.Freq) ~ Bee, BigDataSheet, mean)
BDSRanked$Trial = str_extract(BDSRanked$Bee, ".+?(?=_)")
BDSRanked$Col = str_extract(BDSRanked$Bee, ".+?(?=#)")
BDSRanked = BDSRanked %>% arrange(Col, AntPresence) %>%
group_by(Col) %>%
mutate(Rank = rank(AntPresence))
ggplot(BDSRanked, aes(x = Rank, y = AntPresence)) +
geom_line(data = subset(BDSRanked, QR==1), aes(group = Col, color=Col)) +
geom_point(data = subset(BDSRanked, Queen==1),color="purple",size=3)
Top10 = BDSMeans %>%
group_by(Col) %>%
top_n(10, Degree)
Bottom10 = BDSMeans %>%
group_by(Col) %>%
top_n(10, Degree)
ggplot(BDSMeans, aes(x = QR, y = Degree)) +
geom_line(aes(group = Trial),color="darkgray") +
geom_point(aes(color = Trial), size = 5) +
scale_color_manual(values = wes_palette("Cavalcanti1")) +
xlab("") +
labs(color="Source Colony") +
ylab("Degree") + # Adjust axis labels
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
# legend.position = "none",
panel.grid.major.x = element_line(color = "grey", linetype = "dashed"), # Keep vertical grid lines
panel.grid.minor.x = element_line(color = "grey", linetype = "dotted"), # Keep vertical grid lines
panel.grid.major.y = element_blank(), # Remove horizontal grid lines
panel.grid.minor.y = element_blank(), # Remove horizontal grid lines
axis.line.x = element_line(color = "black", linewidth = 0.5), # Add x-axis line
axis.line.y = element_line(color = "black", linewidth = 0.5), # Add y-axis line
strip.text = element_text(size = 14, face = "bold"), # Make facet plot titles larger and bold
axis.text.y.right = element_blank(), # Remove right y-axis text
axis.ticks.y.right = element_blank(), # Remove right y-axis ticks
aspect.ratio = 1
) +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5))
###AnalyzingTheSheet: OvaryIndex
BDSMeans$OvaryIndex = BDSMeans$AverageOvaryWidth/BDSMeans$AverageWingLength
ggplot(BDSMeans, aes(x=OvaryIndex, y=boutBetween)) +
geom_point() +
geom_smooth(method = lm,
se = T,
linewidth = 1.5,
linetype = 1,
alpha = .7,
aes(color = factor(QR)))+
theme_classic() +
labs(title="Avg Length vs Avg Width")+
scale_color_manual(name ="Bee",
labels = c("QL", "QR"),
values = c("#161414", "#629CC0"))+
facet_wrap(~ QR, nrow = 2) #Made this before, but it's a good sanity check
## `geom_smooth()` using formula = 'y ~ x'
